/*
  Ported to JavaScript&canvas from the paper "Real-Time Fluid Dynamics
  for Games" by Jos Stam(jstam@aw.sgi.com) at GDC2003.
  You can get the original codes and the paper from
  http://www.dgp.toronto.edu/people/stam/reality/Research/pub.html
*/


var N;
var dt;
var diff;
var visc;
var force;
var source;
var dvel = 0;
var u;
var v;
var u_prev;
var v_prev;
var dens;
var dens_prev;
var win_x = 256, win_y = 256;
var mouse_down = new Array(3);
var omx = 0, omy = 0,mx = 0, my = 0;
var canvas;
var ctx;

function array2d(s)
{
  var a = new Array(s);
  for (var i = 0; i < s; i++)
    a[i] = new Array(s);
  return a;
}

function calc_color(c)
{
  var v = (c > 1 ? 255 : Math.floor(c*255)).toColorPart();
  return '#' + v + v + v;
}

function cheap_glquads(bx0,by0,c0,bx1,by1,c1,bx2,by2,c2,bx3,by3,c3) {
  var x = [bx0,bx1,bx2,bx3].map(function(v){return v*win_x});
  var y = [by0,by1,by2,by3].map(function(v){return v*win_y});
  var c = [c0,c1,c2,c3];
  var tx = (x[0]+x[1]+x[2]+x[3]) / 4, ty = (y[0]+y[1]+y[2]+y[3])/4;
  var ac = (c[0] + c[1]+ c[2] + c[3]) / 4;
  var average_color = calc_color(ac);
  var g;
  for (var i = 0; i < 4; i++) {
    g = ctx.createLinearGradient(x[i], y[i], tx, ty);
    g.addColorStop(0,calc_color(c[i]));
    g.addColorStop(1,average_color);
    ctx.fillStyle = g;
    ctx.fillRect( x[i], y[i], tx-x[i], ty-y[i]);
  }
}


function add_source (N,x,s,dt)
{
  for (var i=0 ; i<N+2 ; i++ ) 
    for (var j=0 ; j<N+2 ; j++ )
      x[i][j] += dt*s[i][j];
}

function set_bnd (N,b,x)
{
  var i;
  for ( i=1 ; i<=N ; i++ ) {
    x[0][i] = b==1 ? -x[1][i] : x[1][i];
    x[N+1][i] = b==1 ? -x[N][i] : x[N][i];
    x[i][0] = b==2 ? -x[i][1] : x[i][1];
    x[i][N+1] = b==2 ? -x[i][N] : x[i][N];
  }
  x[0][0] = 0.5*(x[1][0]+x[0][1]);
  x[0][N+1] = 0.5*(x[1][N+1]+x[0][N]);
  x[N+1][0] = 0.5*(x[N][0]+x[N+1][1]);
  x[N+1][N+1] = 0.5*(x[N][N+1]+x[N+1][N]);
}

function lin_solve (N,b,x,x0,a,c)
{
  var i, j, k;
  for ( k=0 ; k<20 ; k++ ) {
    for ( i=1 ; i<=N ; i++ ) {
      for ( j=1 ; j<=N ; j++ ) {
	x[i][j] = (x0[i][j] + a*(x[i-1][j]+x[i+1][j]+x[i][j-1]+x[i][j+1]))/c;
      }
    }
    set_bnd ( N, b, x );
  }
}

function diffuse (N,b,x,x0,diff,dt)
{
  var a=dt*diff*N*N;
  lin_solve ( N, b, x, x0, a, 1+4*a );
}

function advect (N,b,d,d0,u,v,dt)
{
  var i, j, i0, j0, i1, j1;
  var x, y, s0, t0, s1, t1, dt0;

  dt0 = dt*N;
  for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
      x = i-dt0*u[i][j]; y = j-dt0*v[i][j];
      if (x<0.5) x=0.5; if (x>N+0.5) x=N+0.5; i0=Math.floor(x); i1=i0+1;
      if (y<0.5) y=0.5; if (y>N+0.5) y=N+0.5; j0=Math.floor(y); j1=j0+1;
      s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1;
      d[i][j] = s0*(t0*d0[i0][j0]+t1*d0[i0][j1])+s1*(t0*d0[i1][j0]+t1*d0[i1][j1]);
    }
  }
  set_bnd ( N, b, d );
}

function project (N,u,v,p,div)
{
  var i, j;
  for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
      div[i][j] = -0.5*(u[i+1][j]-u[i-1][j]+v[i][j+1]-v[i][j-1])/N;
      p[i][j] = 0;
    }
  }
  set_bnd ( N, 0, div ); set_bnd ( N, 0, p );
  
  lin_solve ( N, 0, p, div, 1, 4 );
  
  for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
      u[i][j] -= 0.5*N*(p[i+1][j]-p[i-1][j]);
      v[i][j] -= 0.5*N*(p[i][j+1]-p[i][j-1]);
    }
  }
  set_bnd ( N, 1, u ); set_bnd ( N, 2, v );
}

function dens_step (N,x,x0,u,v,diff,dt )
{
  add_source ( N, x, x0, dt );
  //SWAP( x0, x ); 
  diffuse ( N, 0, x0, x, diff, dt );
  //SWAP( x0, x );
  advect ( N, 0, x, x0, u, v, dt );
}

function vel_step (N,u,v,u0,v0,visc,dt)
{
  add_source ( N, u, u0, dt ); add_source ( N, v, v0, dt );
  //SWAP ( u0, u );
  diffuse ( N, 1, u0, u, visc, dt );
  //SWAP ( v0, v );
  diffuse ( N, 2, v0, v, visc, dt );
  project ( N, u0, v0, u, v );
  //SWAP ( u0, u ); SWAP ( v0, v );
  advect ( N, 1, u, u0, u0, v0, dt );
  advect ( N, 2, v, v0, u0, v0, dt );
  project ( N, u, v, u0, v0 );
}

function get_from_UI(d,u,v)
{
  var i, j;
  for (i=0 ; i<N+2 ; i++ ) 
    for (j=0 ; j<N+2 ; j++ )
      u[i][j] = v[i][j] = d[i][j] = 0.0;
  
  if (!mouse_down[0] && !mouse_down[1]) return;
  i = Math.floor(( mx /win_x)*N+1);
  j = Math.floor(( my /win_y)*N+1);
  
  if ( i<1 || i>N || j<1 || j>N ) return;
  
  if ( mouse_down[0] ) {
    u[i][j] = force * (mx-omx);
    v[i][j] = force * (omy-my);
  }
  
  if ( mouse_down[1] ) {
    d[i][j] = source;
  }

  omx = mx;
  omy = my;
}

function draw_velocity() {
  var i, j;
  var x, y, h;
  h = 1.0/N;
  ctx.beginPath();
  for ( i=1 ; i<=N ; i++ ) {
    x = (i-0.5)*h;
    for ( j=1 ; j<=N ; j++ ) {
      y = (j-0.5)*h;
      ctx.moveTo( x * win_x, y * win_y);
      ctx.lineTo( (x+u[i][j]) * win_x, (y+v[i][j]) * win_y);
    }
  }
  ctx.stroke();
}

function draw_density()
{
  var i, j;
  var x, y, h, d00, d01, d10, d11;
  h = 1.0/N;
  for ( i=0 ; i<=N ; i++ ) {
    x = (i-0.5)*h;
    for ( j=0 ; j<=N ; j++ ) {
      y = (j-0.5)*h;
      d00 = dens[i][j];
      d01 = dens[i][j+1];
      d10 = dens[i+1][j];
      d11 = dens[i+1][j+1];
      cheap_glquads(x,y,d00,
		    x+h,y,d10,
		    x+h,y+h,d11,
		    x,y+h,d01);
    }
  }
}


function idle() {
  get_from_UI(dens_prev, u_prev, v_prev);
  vel_step(N,u,v,u_prev,v_prev,visc,dt);
  dens_step(N,dens,dens_prev,u,v,diff,dt);
}

function display_func() {
  ctx.clearRect(0,0,win_x,win_y);
  if ( dvel ) draw_velocity ();
  else	      draw_density ();
}

function main() {
  var q = location.search.toQueryParams();
  N = q.N ? parseInt(q.N) : 16;
  dt = q.dt ? parseFloat(q.dt) : 0.1;
  diff = q.diff ? parseFloat(q.diff) : 0.0;
  visc = q.visc ? parseFloat(q.visc) : 0.0;
  force = q.force ? parseFloat(q.force) : 5.0;
  source = q.source ? parseFloat(q.source) : 100.0;
  u = array2d(N+2);
  v = array2d(N+2);
  u_prev = array2d(N+2);
  v_prev = array2d(N+2);
  dens = array2d(N+2);
  dens_prev = array2d(N+2);

  canvas = $('main');
  canvas.onmousemove = motion_func;
  canvas.onmousedown = mouse_func;
  canvas.onmouseup = mouse_func2;
  ctx = canvas.getContext('2d');

  for (var i=0 ; i<N+2 ; i++ ) 
    for (var j=0 ; j<N+2 ; j++ )
      u[i][j] = v[i][j] = u_prev[i][j] = v_prev[i][j] = dens[i][j] = dens_prev[i][j] = 0.0;  
  setInterval(display_func,50);
  setInterval(idle,50);
}

function mouse_func(e) {
  if (!e) var e = window.event;
  omx = mx = e.clientX;
  omx = my = e.clientY;
  if(e.shiftKey)
    mouse_down[0] = true;
  else 
    mouse_down[1] = true;
}

function mouse_func2(e) {
  if (!e) var e = window.event;
  if(e.shiftKey)
    mouse_down[0] = false;
  else
    mouse_down[1] = false;
}

function motion_func(e) {
  if (!e) var e = window.event;
  mx = e.clientX;
  my = e.clientY;
}

window.onload = main;

